/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  
  static int solve_ldl_nopiv_1_flat( int fs, int ss,
                                     linal::Flat_ ATL, linal::Flat_ ABL, 
                                     linal::Flat_ t,   linal::Flat_ b );

  static int solve_ldl_nopiv_1_hier( int fs, int ss,
                                     linal::Hier_ ATL, linal::Hier_ ABL, 
                                     linal::Hier_ t,   linal::Hier_ b );
  
  static int solve_ldl_nopiv_2_flat( int fs, int ss,
                                     linal::Flat_ ATL, linal::Flat_ ATR,
                                     linal::Flat_ ABL,
                                     linal::Flat_ t,   linal::Flat_ b );

  static int solve_ldl_nopiv_2_hier( int fs, int ss,
                                     linal::Hier_ ATL, linal::Hier_ ATR,
                                     linal::Hier_ ABL,
                                     // ***************
                                     linal::Flat_ FTL, linal::Flat_ FTR,
                                     linal::Flat_ FBL,
                                     // ***************
                                     linal::Hier_ t,   linal::Hier_ b );
    
  void Matrix_FLA_::solve_ldl_nopiv_1_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_1_hier( this->fs, this->ss,
                            this->hier.ATL, this->hier.ABL,
                            this->hier.xt,  this->hier.xb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_1_flat( this->fs, this->ss,
                            this->flat.ATL, this->flat.ABL,
                            this->flat.xt,  this->flat.xb );

#endif
  }

  void Matrix_FLA_::solve_ldl_nopiv_2_x() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_ldl_nopiv_2_hier( this->fs, this->ss,
                            this->hier.ATL, this->hier.ATR,
                            this->hier.ABL,
                            // ***************                                   
                            this->flat.ATL, this->flat.ATR,
                            this->flat.ABL,
                            // ***************                                   
                            this->hier.xt,  this->hier.xb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_2_flat( this->fs, this->ss,
                            this->flat.ATL, this->flat.ATR,
                            this->flat.ABL,
                            this->flat.xt,  this->flat.xb );
#endif
  }    

  void Matrix_FLA_::solve_ldl_nopiv_1_r() {
#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_1_hier( this->fs, this->ss,
                            this->hier.ATL, this->hier.ABL,
                            this->hier.rt,  this->hier.rb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_1_flat( this->fs, this->ss,
                            this->flat.ATL, this->flat.ABL,
                            this->flat.rt,  this->flat.rb );
#endif
  }

  void Matrix_FLA_::solve_ldl_nopiv_2_r() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    solve_ldl_nopiv_2_hier( this->fs, this->ss,
                            this->hier.ATL, this->hier.ATR,
                            this->hier.ABL,
                            // ***************                                   
                            this->flat.ATL, this->flat.ATR,
                            this->flat.ABL,
                            // ***************                                   
                            this->hier.rt,  this->hier.rb );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix 
    // ----------------------------------------------------------
    solve_ldl_nopiv_2_flat( this->fs, this->ss,
                            this->flat.ATL, this->flat.ATR,
                            this->flat.ABL,
                            this->flat.rt,  this->flat.rb );
#endif

  }    
  					 
  static inline int solve_ldl_nopiv_1_flat( int fs, int ss,
                                            linal::Flat_ ATL, linal::Flat_ ABL, 
                                            linal::Flat_ t,   linal::Flat_ b ) {
    if (fs) 
      FLA_Trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
                FLA_NO_TRANSPOSE,
                FLA_UNIT_DIAG, FLA_ONE,
                ~ATL, ~t );
    
    if (fs && ss)
      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ABL, ~t,
                FLA_ONE, ~b );
    
    return true;
  }

  static inline int solve_ldl_nopiv_1_hier( int fs, int ss,
                                            linal::Hier_ ATL, linal::Hier_ ABL, 
                                            linal::Hier_ t,   linal::Hier_ b ) {

    if (fs) 
      linal::trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
                          FLA_NO_TRANSPOSE,
                          FLA_UNIT_DIAG, FLA_ONE,
                          ATL, t );

    if (fs && ss) 
      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                          FLA_MINUS_ONE, ABL, t,
                          FLA_ONE, b);

    return true;
  }

  static inline int solve_ldl_nopiv_2_flat( int fs, int ss,
                                            linal::Flat_ ATL, linal::Flat_ ATR,
                                            linal::Flat_ ABL, 
                                            linal::Flat_ t,   linal::Flat_ b ) {
    if (ss && fs) { 

      ATR.create( ABL.get_datatype(), ABL.get_n(), ABL.get_m() );
      
      FLA_Copyt( FLA_CONJ_TRANSPOSE, ~ABL, ~ATR );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE,
                                              ~ATL, ~ATR );

      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ATR, ~b,
                FLA_ONE, ~t );

      ATR.free();
    }

    if (fs) {
      FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~ATL );
   
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              FLA_LOWER_TRIANGULAR, ~ATL );

      FLA_Trsm( FLA_LEFT, FLA_UPPER_TRIANGULAR,
                FLA_NO_TRANSPOSE,
                FLA_NONUNIT_DIAG, FLA_ONE,
                ~ATL, ~t );
    }
    return true;
  }

  static int solve_ldl_nopiv_2_hier( int fs, int ss,
                                     linal::Hier_ ATL, linal::Hier_ ATR, 
                                     linal::Hier_ ABL, 
                                     // ***************
                                     linal::Flat_ FTL, linal::Flat_ FTR, 
                                     linal::Flat_ FBL, 
                                     // ***************
                                     linal::Hier_ t,   linal::Hier_ b ) {
    if (ss && fs) {
      int bmn = get_hier_blocksize();

      FTR.create( FBL.get_datatype(), FBL.get_n(), FBL.get_m() );
      ATR.create( FTR, bmn, bmn );
      
      FLA_Copyt( FLA_TRANSPOSE, ~FBL, ~FTR );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              ~FTL, ~FTR );

      linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                          FLA_MINUS_ONE, ATR, b,
                          FLA_ONE, t );

      ATR.free();
      FTR.free();
    }

    if (fs) {
      FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~FTL );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              FLA_UPPER_TRIANGULAR, ~FTL );

      linal::trsm( FLA_LEFT, FLA_UPPER_TRIANGULAR,
                          FLA_NO_TRANSPOSE,
                          FLA_NONUNIT_DIAG, FLA_ONE,
                          ATL, t );
    }

    return true;
  }
}
